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A  modification  of  a  previously  developed  computer  model  of  snow- 
melt  provides  for  modeling  intermittent  snowpacks,  and  is  believed  to  be 
a  more  generalized  model  than  the  original  program.  The  modified 
program  SNOWMELT  is  dependent  on  four  daily  input  variables, 
maximum  and  minimum  temperatures,  precipitation,  and  shortwave 
radiation  or  percent  cloud  cover.  Initializing  the  model  requires  limited 
knowledge  of  local  watershed  and  snowpack  parameters.  Model  verifica- 
tion on  seven  experimental  watersheds  in  Arizona  proved  satisfactory. 

Keywords:  Computer  models,  intermittent  snowpacks,  simulation  analy- 
sis, snowmelt. 
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Computer  Simulation  of  Snowmelt 


Rhey  M.  Solomon,  Peter  F.  Ffolliott,  Malchus  B.  Baker,  Jr.,  and  J.  R.  Thompson 


Introduction 

Researchers  have  studied  snow  in  various  ways  to 
more  fully  understand  processes  involved  in  snow 
accumulation,  melt,  and  its  resultant  runoff.  To 
investigate  interrelationships  of  the  many  variables 
affecting  snowmelt  would  require  large  amounts  of 
data,  many  parameters,  and  numerous  calculations. 
To  circumvent  many  of  these  problems,  computer 
models  have  been  synthesized  that  are  designed  to 
simulate  the  melting  of  snowpacks  (Leaf  and  Brink 
1973a,  1973b;  U.S.  Army,  Corps  of  Engineers  1972). 

For  computer  models  to  be  of  practical  use,  the 
input  variables  required  for  the  models  should  be 
easily  obtained,  and  the  models  must  have  utility  for 
the  land  manager. 

Leaf  and  Brink  (1973a)  developed  a  model  (MELT- 
MOD)  intended  to  be  used  with  input  variables 
sensitive  to  vegetation  manipulations.  MELTMOD 
was  developed  for  Colorado  subalpine  forests,  with 
verification  data  collected  from  the  667-acre  Dead- 
horse  Creek  watershed  on  the  Fraser  Experimental 
Forest  (Leaf  1971,  Leaf  and  Brink  1972).  Thus,  the 
model  was  limited  in  its  scope  of  possible  melt 
regimes,  and  by  snowpack  conditions  that  might  be 
encountered  in  other  areas  of  the  country. 

MELTMOD  was  evaluated  as  a  model  for  estimat- 
ing daily  and  weekly  snowmelt  rates  for  snowpack 
conditions  found  in  the  Southwest,  and  more  specif- 
ically as  a  tool  for  characterizing  runoff  efficiencies  of 
small  Arizona  watersheds  (Solomon  et  al.  1975). 
MELTMOD  was  found  to  be  unsatisfactory,  however, 
principally  because  of  three  shortcomings:  the  in- 
ability of  the  model  to  respond  quickly  to  alternate 
freezing  and  melting  patterns  during  the  accumula- 
tion and  melt  phases;  the  model  tended  to  keep  the 
snowpack  excessively  "cold"  during  the  accumulation 
phase;  and  some  of  the  input  variables  necessary  for 
the  model,  daily  solar  radiation  in  particular,  were 
not  always  available.  Because  of  these  inadequacies  in 
modeling  intermittent  snowpacks,  MELTMOD  was 
modified  to  account  for  intermittent  snowpack  condi- 
tions which  are  often  found  in  the  Southwest.  The 
intent  of  this  modification  was  to  produce  a  general 
snowmelt  model,  applicable  to  simulating  both  con- 
tinuous and  intermittent  snowpacks.  However,  all 
tests  of  the  modified  model  to  date  have  pertained  to 
intermittent  snowpacks,  and  no  attempt  has  yet  been 
made  to  verify  the  new  model's  applicability  for 
simulating  continuous  snowpacks. 


Program  Description 

Basically,  the  modified  computer  model  SNOW- 
MELT simulates:  winter  snow  accumulation;  energy 
balance;  snowpack  conditions;  and  resultant  snow- 
melt, with  a  sensitivity  to  slope,  aspect,  temperature, 
shortwave  radiation,  and  forest  cover  density. 
SNOWMELT  simulates  the  same  basic  processes  that 
were  common  to  MELTMOD:  the  segregation  of 
precipitation  into  the  snow  and  rain  components;  the 
melting  process,  including  the  energy  budget;  and 
snowpack  conditions  in  terms  of  energy  level  and  free 
water  requirements.  Modifications  of  MELTMOD 
were  necessary,  however,  to  better  simulate  the  inter- 
mittent snowpacks  common  in  Arizona.  Many  of  the 
subroutines  in  the  MELTMOD  model  have  not  been 
altered,  and  their  flow  paths  will  not  be  discussed; 
subroutines  that  have  been  eliminated,  added,  or 
modified  will  be  reviewed.3 

Five  alterations  of  MELTMOD  were  necessary: 
segregation  of  rain  and  snow  inputs,  synthesis  of  solar 
radiation  data,  alteration  of  snowpack  radiating 
temperature,  design  of  a  thermal  diffusion  technique, 
and  an  accounting  for  areal  extent  during  snowpack 
ablation. 

The  modified  program  SNOWMELT  is  dependent 
on  only  four  daily  input  variables:  (1)  maximum 
temperature,  (2)  minimum  temperature,  (3)  precipi- 
tation, and  (4)  shortwave  radiation  or  percent  cloud 
cover.  In  addition,  the  following  initializing  values 
are  required:  (1)  initial  snowpack  temperature,  ^2)  a 
solar  radiation  transmissivity  coefficient,  (3)  forest 
cover  density,  (4)  initial  snowpack  water  equivalent, 
(5)  a  threshold  temperature  value  for  use  in  calculat- 
ing snowpack  reflectivity,  (6)  mean  slope  and  aspect 
of  the  watershed,  (7)  latitude,  (8)  an  atmospheric 
absorption  coefficient,  and  (9)  a  time  interval  for 
hour  angle  changes  in  subroutine  SOLAR.  If  daily 
solar  radiation  values  are  available,  they  are  corrected 
for  slope  and  aspect  prior  to  use  in  SNOWMELT. 
Therefore,  variables  (6),  (7),  (8),  and  (9)  are  not 
necessary  in  SOLAR,  as  they  are  only  used  in  calcu- 
lating a  daily  insolation  incident  upon  the  forest 
canopy. 


All  subroutines  have  been  programed  for  the  CDC  6400  computer 
at  the  University  of  Arizona,  Tucson. 
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Segregating  Precipitation  Into 
Rain  and  Snow  Components 

Separation  of  precipitation  into  snow  and  rain 
components  is  one  of  the  first  steps  in  the  flow  of  the 
model  (fig.  1).  The  original  MELTMOD  program 
used  the  daily  minimum  temperature  as  an  index  in 
segregating  precipitation  into  snow,  rain,  or  mixed 
events.  If  the  daily  minimum  temperature  fell  below 
32°F  or  the  daily  maximum  temperature  was  less 
than  35°F,  any  precipitation  fell  as  snow.  We  altered 
MELTMOD  guidelines  to  consider  precipitation  as 


all  snow  if  the  maximum  temperature  drops  below 
40°  F  and  the  minimum  daily  temperature  never 
exceeds  35° F.  Criteria  for  mixed  snow  and  rain 
events  are  altered  only  slightly.  Precipitation  falls  as 
rain  when  the  minimum  temperature  exceeds  35°F; 
mixed  rain  and  snow  falls  when  the  minimum  tem- 
perature drops  below  35°F,  and  the  maximum  tem- 
perature exceeds  40°F.  The  mixed-event  partitioning 
is  performed  through  the  same  equation  as  used  in 
MELTMOD  (Leaf  and  Brink  1973a).  The  reasons  for 
these  changes  came  from  examination  of  winter- 
precipitation  data  from  Flagstaff,  Arizona  (Beschta 
1974). 


^START^ 


Initialize 

variables 

Read  a  card 


Yes 


Determine  type  of  precip.  and 
its  energy  and  mass  contributions 


Calculate  the  daily  energy  budget 


Determine  snowpack  temperature 
with  the  diffusion  model  (DIFTEM) 


Calculate  effects  of  energy  inputs 
on  the  snowpack 


Store  results 


Print  daily  results  (water  eq.,  temp.,  etc.) 


^END^ 


Figure  1.— General  flow  chart  for  program  SNOWMELT. 
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Synthesis  of  Solar  Radiation  Data 

Lack  of  daily  solar  radiation  data  prompted  the 
incorporation  of  subroutines  SOLAR  and  CLOUD 
into  the  program  (see  appendix). 

Subroutine  SOLAR  computes  potential  daily  solar 
insolation  inputs  for  uniform  increments  of  time, 
routes  the  insolation  through  the  atmosphere,  and 
corrects  for  slope  and  aspect.  Incremental  values  are 
then  summed  to  arrive  at  total  insolation  for  a  day 
(fig.  2). 


Initialize    counters  and  variables 


Compute  the  number  of  days  since  winter  solstice 


Add  an  increment  to  sun's  path 


Calculate  optical  path  from  zenith  angle 


Calculate  direct  beam  radiation  on  a  flat  surface 


Calculate  ratio  of  insolation  that  is  diffuse  beam 
I 

Calculate  direct  and  diffuse  beam  accounting  for  slope 
I 

Add  direct  and  diffuse  beam  components 


Add  the  insolation  for  this  interval  to  the  daily  accumulated  insolation 


T 


Call  CLOUD 


Figure  2.  — Flow  chart  for  subroutine  SOLAR. 


In  this  subroutine,  direct  beam  solar  insolation  is 
obtained  from  equations  presented  by  Frank  and  Lee 
(1966).  These  values  are  then  routed  through  the 
atmosphere  by  an  experimentally  determined  "extinc- 
tion" coefficient.  For  Arizona  conditions,  this  coef- 
ficient was  found  to  be  0.13.  Diffuse-beam  solar 
insolation  is  determined  from  knowledge  of  solar 
altitude  and  direct-beam  solar  insolation  by  equations 
empirically  derived  from  tables  given  by  Reifsnyder 
and  Lull  (1965).  A  cloudless  atmosphere  is  assumed. 


Potential  daily  solar  insolation  is  calculated  every 
fifth  day,  since  changes  over  a  5-day  period  are  small. 

Once  a  method  for  determining  values  for  clear 
days  was  derived,  some  adjustment  was  necessary  for 
cloudy  days  (subroutine  CLOUD).  An  equation  pre- 
sented by  Gates  (1962),  using  mean  cloud  cover  and  a 
coefficient  dependent  on  latitude,  proved  to  be 
satisfactory: 

Q  =  Q0[l  -  (l-K)C] 

where: 

Q0  =  total  energy  received  without  accounting  for 
cloud  cover; 

K  =  a  cloud  cover  parameter  which  is  a  function  of 
latitude,  given  a  value  0.38  at  30°,  0.41  at  40°, 
and  0.42  at  50°;  and 

C  =  average  daily  proportion  of  clouded  to  clear 
skies. 

The  cloud  coefficients  (K)  were  increased  slightly 
from  those  shown  by  Gates  because  the  equation 
tends  to  overestimate  radiation  for  days  with  consid- 
erable cloud  cover  and  underestimate  for  days  with  50 
to  70  percent  cloud  cover.  Testing  calculated  values 
against  3  months  of  Flagstaff  source  data  (January 
through  March  1974)  yielded  an  average  daily  value 
that  was  within  96  percent  of  the  actual  daily  insola- 
tion and  had  a  correlation  coefficient  of  0.87. 


Snowpack  Radiating  Temperature 


When  a  caloric  deficit  existed  in  MELTMOD  (a 
snowpack  below  0°C),  the  radiating  temperature  of 
the  snowpack  was  made  equal  to  the  average  ambient 
air  temperature  or  0°C,  whichever  was  less.  If  there 
was  no  caloric  deficit,  the  radiating  temperature  was 
made  equal  to  the  minimum  air  temperature  or  0°C, 
whichever  was  less. 

Changing  the  radiating  temperature  with  snowpack 
conditions  was  necessary  because  of  Colorado  snow- 
pack conditions.  In  Colorado,  snowpacks  normally 
increase  in  water  content  throughout  the  accumula- 
tion phase,  with  little  or  no  melting  (Leaf  and  Brink 
1973a).  Therefore,  to  keep  the  pack  "cold"  during 
the  accumulation  phase,  thermal  exchanges  were 
controlled  by  a  diffusion  model.  By  using  average 
ambient  air  temperature,  the  snowpack  remained 
below  0°C  until  the  melt  season,  at  which  time  the 
model  could  be  reinitialized  by  making  the  snowpack 
temperature  0°C.  This  is  satisfactory  for  an  area  with 
distinct  accumulation  and  melt  season,  but  where 
periods  of  melt  and  accumulation  alternate  through- 
out the  winter,  as  in  Arizona,  this  approach  is  in- 
adequate. Therefore,  the  radiating  snowpack  temper- 
ature was  altered  to  allow  only  the  minimum  air 
temperature  or  0°C,  whichever  was  less,  to  be  used  as 
the  radiating  snowpack  temperature. 
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Thermal  Diffusion  Subroutine 

MELTMOD  used  a  modified  energy  budget  during 
periods  of  negative  energy  inputs.  The  routines  con- 
trolling the  use  of  the  thermal  diffusion  model, 
principally  DIFMOD,  LINK,  and  RDPACK,  were 
designed  to  control  the  snowpack  temperature  during 
the  winter  months  when  energy  losses  from  the  snow- 
pack  kept  its  temperature  below  0°C.  Additionally, 
the  DIFMOD  routine  was  only  mathematically  stable 
for  snowpacks  greater  than  4.7  inches  in  depth.  As  a 
result,  all  three  subroutines  were  discarded  and  a 
single  subroutine  (DIFTEM)  was  written  for  Arizona 
conditions.  A  flow  chart  of  DIFTEM  is  shown  in 
figure  3,  and  the  subroutine  is  given  in  the  appendix. 
A  central  differencing  technique  (Carnahan  et  al. 
1969)  was  modified  and  adjusted  to  program  SNOW- 
MELT. 

In  the  forward  differencing  procedure  used  by  Leaf 
and  Brink  (1973a)  in  the  development  of  MELT- 
MOD,  the  value  of  U:  :  is  dependent  upon  \ii  i  ■ 
i»  J  J"i 

Uj        and  u-_|_j  j  j  (fig.  4).  Only  those  values  of  u 

within  the  pyramidal  area  A  can  have  any  influence 
on  the  value  of  U:  •,  whereas  it  is  known  that  points 

in  time  earlier  than  tj  in  area  B  also  influence  the 

value  of  U:       In  the  implicit  central  differencing 
*>  J 

method  we  use,  values  of  u-  •  are  evaluated  at  the 

*>  J 

advanced  point  in  time  t j  _|_  j  instead  of  tj ,  as  in  the 


^  DIFTEM^ 


Normalize   minimum  air  temp  and  0°C  to  a  0  to  1.0  scale 


No 


Calculate   snowpack  density 


Determine   an  average  snowpack  temp. 


Normalize   temperatures  throughout  the  snowpack 


Calculate  new  temperature  profiles   over  succesive  time-steps 
throughout   a  series  of  numerical  approximations 


Establish  the  array  of  final  pack  temperatures   to  be  used 
as  the  array  and  calorie  deficit  for  the  next  day 


Determine  the  calorie  deficit  and  average  pack  temp. 


1_ 


Figure  3.  — Flow  chart  for  subroutine  DIFTEM. 

forward  differencing  procedure.  This  provides  a  wider 
base  from  which  u5  :  can  be  estimated. 
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Snowpack  Conditions 

Because  snowpack  ablation  and  areal  depletion  are 
not  uniform  throughout  a  forest,  snow  cover  becomes 
discontinuous  toward  the  end  of  the  melt  season 
(Brown  and  Dunford  1956,  Leaf  1971).  The  original 
MELTMOD  model  did  not  account  for  this  phenom- 
enon, and  as  a  result  synthesized  snowpack  ablation 
rates  were  often  excessive  near  the  end  of  the  melt 
season;  100  percent  ablation  was  indicated  before  the 
end  of  generated  runoff. 

Snowpack  water-equivalent  measurements  taken 
on  small  experimental  watersheds  in  Arizona  were 
utilized  in  determining  snowpack  cover  percent  as  a 
function  of  mean  watershed  snowpack  water-equiva- 
lent. For  each  sampling  date,  the  number  of  sample 
points  with  no  snow  were  ratioed  to  the  total  points 
surveyed,  and  the  mean  snowpack  water-equivalents 
were  calculated  for  that  date.  Graphical  analysis  of 
the  data  revealed  that  only  snowpack  water-equiva- 
lents under  3  inches  showed  deviation  from  100  per- 
cent snow  cover.  Through  regression  analysis,  27 
snowpack  water-equivalent  measurements  were  re- 
lated to  their  corresponding  percentages  of  snow 
cover.  The  correlation  coefficient  for  the  following 
equation  was  0.93: 

percent  cover  =  31.9  (snowpack  water-equivalent) 

This  equation  is  only  used  for  snowpacks  contain- 
ing 3  inches  or  less  of  water.  Therefore,  when  the 
mean  watershed  snowpack  water-equivalent  is  below 
3  inches,  the  calculated  melt  is  reduced  by  the  per- 
centage of  area  not  covered  by  snow.  This  technique 
is  one  method  of  simulating  observations  made  by 
other  researchers  (Gary  and  Coltharp  1967,  Ffolliott 
and  Hansen  1968,  Leaf  1971). 


Model  Verification 

The  model  was  verified  by  processing  source  data 
collected  from  seven  experimental  watersheds  ranging 
from  20  to  467  acres.  These  watersheds  were  located 
across  Arizona  to  obtain  an  array  of  vegetative, 
physiographic,  and  climatic  conditions  found  in  the 
Southwest. 


Watershed  Descriptions 

Beaver  Creek  watershed  15,  an  area  of  163  acres,  is 
located  approximatey  20  miles  south  of  Flagstaff.  The 
forest  overstory  is  predominantly  ponderosa  pine, 
with  an  intermixing  of  Gambel  oak  and  alligator 
juniper.  Soils  are  derived  from  volcanic  parent 
materials;  elevations  range  from  6,735  to  7,160  feet. 
Annual  precipitation  averages  28  inches,  one-third  of 
which  comes  during  November  through  April.  Princi- 
pal, drainage  is  toward  the  southeast. 


Beaver  Creek  watershed  17,  about  299  acres  in 
size,  is  2  miles  southeast  of  watershed  15.  Ponderosa 
pine  dominates  the  forest  overstory,  although  Gambel 
oak  and  alligator  juniper  occur  in  scattered  clumps. 
Soils  are  derived  from  volcanics;  elevations  range 
from  6,830  to  7,200  feet.  Principal  drainage  is  toward 
the  southwest. 

Heber  watersheds  1  and  2,  approximately  20  and 
28  acres,  respectively,  are  located  13  miles  southeast 
of  Heber,  Arizona.  The  primary  forest  overstory  of 
both  watersheds  is  ponderosa  pine,  with  an  intermix- 
ing of  Gambel  oak,  Douglas-fir,  and  white  fir.  Soils 
are  derived  from  alluvial  parent  materials;  elevations 
range  from  7,400  to  7,700  feet.  Annual  precipitation 
is  approximately  22  inches,  half  of  which  comes 
during  the  winter.  Both  watersheds  have  principal 
drainage  to  the  south. 

Heber  watersheds  3  and  4,  60  and  61  acres,  respec- 
tively, are  located  6  miles  south  of  Heber,  Arizona. 
Principal  forest  overstory  species  are  ponderosa  pine 
and  Gambel  oak,  with  an  intermixing  of  alligator 
juniper.  Soils  are  derived  from  sandstone  parent 
materials;  elevations  range  from  6,900  to  7,050  feet. 
Both  watersheds  drain  toward  the  north,  with  annual 
precipitation  patterns  and  amounts  similar  to  the 
Heber  watersheds  1  and  2. 

The  North  Fork  of  Thomas  Creek,  a  467-acre 
watershed,  is  located  20  miles  south  of  Alpine.  The 
mixed  conifer  forest  overstory  is  predominantly 
Douglas-fir,  white  fir,  corkbark  fir,  Engelmann 
spruce,  and  quaking  aspen  on  north  slopes,  with 
ponderosa  pine  and  Gambel  oak  on  south  slopes. 
Soils  are  derived  from  basalt  parent  materials; 
elevations  range  from  8,400  to  9,150  feet.  One  pri- 
mary stream  channel  traverses  the  watershed  from 
southwest  to  northeast.  Annual  precipitation  averages 
27  inches,  one-third  of  which  comes  during  November 
through  April. 

Assessment  of  SNOWMELT 


Two  input  variables  used  in  SNOWMELT — trans- 
missivity  and  forest  cover  density — were  found  to 
have  a  marked  influence  on  net  daily  energy  budgets. 
An  alteration  of  as  little  as  1  percent  in  transmissivity 
can  produce  a  change  of  0.2^inch  in  predicted  snow- 
pack water  equivalent  by  the  end  of  the  snowmelt 
season.  Values  of  forest  cover  density  also  affect  melt 
rates,  but  not  to  the  same  extent  as  transmissivity. 
Forest  cover  density  also  plays  an  important  role  in 
maintaining  a  "cold"  snowpack  during  the  seasonal 
accumulation  phase.  While  transmissivity  and  forest 
cover  density  are  the  most  important  variables  in 
simulating  the  melting  of  snowpacks,  it  should  be 
noted  that  these  are  also  the  most  difficult  variables 
to  estimate. 

The  use  of  the  new  subroutine  for  periods  of  nega- 
tive energy  links  proved  to  be  satisfactory  in  not 
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permitting  the  snowpack  temperature  to  fall  below 
-5°C  for  more  than  a  1-  or  2-day  period.  It  also 
allowed  quick  recovery  of  the  snowpack  to  isothermal 
conditions  during  periods  of  positive  energy  inputs. 

Difficulties  were  encountered  in  the  use  of  extrapo- 
lated temperature,  cloud  cover,  and  precipitation 
data  due  to  failure  of  equipment  located  near  the 
watersheds.  In  some  cases,  data  were  extrapolated 
over  50  miles.  During  some  of  these  periods,  separa- 
tion of  precipitation  into  snow  and  rain  components 
was  believed  to  have  been  inaccurate. 

It  is  difficult  to  statistically  assess  the  "goodness" 
of  any  snowmelt  prediction  in  time-series  models 
because  of  serial  dependence  within  the  model  (Chow 
1964,  Yevjevich  1972).  Therefore,  of  the  seasons 
modeled  (13  seasons  of  record  on  seven  watersheds), 
three  synthesized  snowmelt  seasons  are  graphically 
presented  in  figure  5:  (A)  a  season  of  "good"  fit 
(watershed  17,  1969);  (B)  a  season  of  "average"  fit 
(watershed  15,  1968);  and  (C)  a  season  depicting  a 
"poor"  fit  (North  Fork  of  Thomas  Creek,  1973), 
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Figure  5.— Three  seasons  of  snowmelt  synthesized  from 
program  SNOWMELT. 


Other  peculiarities  of  the  model  were  also  noted.  In 
the  initial  trials  the  North  Fork  of  Thomas  Creek  was 
modeled  as  an  entity,  but  program  SNOWMELT 
failed  to  produce  melt  patterns  corresponding  to 
measured  snowmelt   water-equivalents.  Therefore, 


because  of  its  steep  north-  and  south-facing  slopes, 
the  North  Fork  of  Thomas  Creek  was  divided  into  two 
contributing  subwatersheds.  Each  area  was  then 
assessed  as  if  it  were  a  separate  watershed.  By  this 
technique,  melt  regimes  were  generated  that  passed 
"near"  points  of  observed  snowmelt  water-equivalent 
for  each  slope.  It  was  concluded  that  watersheds 
having  steep  slopes  facing  north  and  south  should  be 
modeled  as  subwatersheds,  with  each  slope  assessed 
as  a  separate  hydrologic  unit.  No  modeling  difficul- 
ties were  encountered  when  watersheds  containing 
steep  slopes  facing  east-west  (Heber  watersheds  1  and 
2)  were  modeled  as  one  contributing  unit. 


Conclusions 


SNOWMELT  is  believed  to  be  a  more  generalized 
model  than  the  original  MELTMOD  program,  and  it 
appears  to  have  satisfactorily  simulated  intermittent 
snowpack  conditions  and  snowmelt  regimes.  Initializ- 
ing of  the  model  requires  only  limited  knowledge  of 
local  watershed  and  snowpack  parameters.  The  only 
driving  variables  required  for  SNOWMELT  are  daily 
values  of  maximum  and  minimum  air  temperature, 
precipitation,  and  cloud  cover  or  solar  radiation. 
With  some  additional  refinements  (particularly  in 
snow  density  predictions)  and  as  data  become  avail- 
able to  more  accurately  characterize  snow  cover 
densities  and  transmissivities,  SNOWMELT  can  be 
used  over  a  wide  range  of  snowpack  regimes. 
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APPENDIX 


Subroutine  SOLAR 


THIS  SUBROUTINE  IS  DESIGNEt,  TO  COMPUTE  THE  AMOUNT  OF  INCOMING 
SOLAR  RADIATION  INCIDENT  ON  A  GIVEN  SURFACE  DURING  CLEAR  SKY. 
DICTIONARY 

ABSK   -  ABSORPTION   COEFFICIENT   FOR   THE  ATMOSPHERE 
ASPECT  -  AZIMUTH  ANGLE,    DEGREES      FROM  NORTH 
COSZ  -  COSINE   OF   ZENITH   ANGLE ( Z ) 
D  -   INCOMING   DIFFUSE   EEAM   RADIATION,  LYS/MIN 
DELH  -  CHANGE    IN  HOUR  AN GL E ( H ) , M  INUT ES 

DRATIO  -  RATIO  OF  DIFFUSE  TO  DIRECT  BEAM  RA0IATI0N  AS  A  FUNCTION 
OF  Z 

H   -  HOUR  ANGLE 

HI    -  HOURS 

H2    -  MINUTES 

OP  -  OPTICAL  PATH 

0  -  INCOMING  DIRECT   BEAM   RADIATION,  LYS/MIN 

0P0T  -   POTENTIAL   DIRECT   BEAM   RADIATION   ON    A   NORMAL   SURFACE  ABOVE 

THE  E ARTH "S   ATMOSPHERE,  LYS/MIN 

S    -  SOLAR    CONSTANT,    2  LYS/MIN 

SLOPE  -  SLOPE,  PERCENT 

XDEC  -  SOLAR  DECLINATION 

XLAT  -  LATITUDE,  DEGREES 

Z   -  ZENITH  ANGLE 

COMMON   ACTDATE, ACTUAL(21 ) , A VETE M C , A  V ET EMF , B AS T E MF . E NGB AL . HOL DC  A P 
COMMON  CALAIR, CALDEF, CALORIE, CALSN0W,C0VDEN, DATE, DATES(3), DEN 
COMMON   DENSITY,IFIRST,F00TN0T(16),FREEWAT.IDATE(3  72),ISN0W,ITABLE 
COMMON  K0UNT»LASTUSD, LINES, OBSWEQV»RADIN,RADLWN,RADSWN, REFLECT 
COMMON  PACKTEM,PARTICE,PASTINT,PLOTOBS,PLOTWE,FRECIP,PREWEQV 
COMMON  SN0MELT,S0BSEQV(372),SPRECIP(3  72),SPREQV(3  72),SUBTITL(8) 
COMMON  TOTPREC,TITLE(b),THRSHLD,TEMPMIN,TEMPMAX,TCOEFF,USEMEAN 
COMMON  XMAX,  USE  POT, SLOPE,  ASPECT, XLAT, DELH, ABSK, CLC0V,ZAT,SUM<, 
COMMON  KKPP,XLATE 

INTEGER   ACTDATE, DATE, DATES, FOOT NOT.PASTI NT, PLOTOBS,PL0TWE,SUBTITL 
INTEGER   TITLE, US E M E AN, U S E POT 
DIMENSION  ZCLAN ( 12 ) 

DATA   ZC LAN/0. 0,31. 0,59. 0,90. 0,120. 0,151. 0.181. 0,212. 0,2*3.0, 
+  273. 0,30<.. 0,33*1.0/ 
PI   •  3.1*159265 
S   -  2.0 

I F ( KKPP.NE. 5)    GO  TO  902 

DELT   »  DELH 

K    ■  (2*.*60.)/DELH 

DELH   *  (0ELH/60.)*15./180.*PI 

XLATE   »  XLAT 

XLAT  •  <XLAT/1B0)*PI 

SLOPE   «  ATANISLOPE/100) 

ASPECT  •-(ASPECT/180. )*PI 


902   KKPP   «   KKPP   +  1 

IFIKKPP.LT. 5)    GO  TO  103 

KKPP    «  0 

H      »  PI 

SUM1  »0. 

SUM2  -0. 

SUM3  *0. 

IDB  »DATES(1> 

X    »  ZCLAN(IDB) 

DAY   *  DATES ( 2 )    +  X 

IF(0AY.GE.355  )   GO  TO  11* 

DAY   =»    DAY   +  11 

GO   TO  115 
11*   DAY   *   DAY   -  355 

115   XOEC   •  -.*1015*(CQS(DAY+PI/1S3>> 

X3   =■   SIN( SLOPE )*CQS (ASPECT) *C0S( XLAT )+C0S ( SLOPE )*SIN( XLAT) 
X*   «  ((l.-X3**2)**O.5)*C0S(X0EC) 

X5    •  ATAN((SIN(SLOPE)*SIN(ASPECT))/(COS(SLOPE)*COS(XLAT)-SIN(SLnPE 
•)*COS(ASPECT)*SIN(XLAT))) 
X6   «  X3*SIN(XDEC) 
DO   100  I-l.K 
H   ■  H— DEL  H 

COSZ   ■  SIN( XLAT ) *S IN1XDEC ) +COS ( XLAT ) *C0S ( XDEC )*COS ( H) 

IF (COSZ. LT. 0. )    GO  TO  100 

Z»AC0S(C0SZ)*(180./PI) 

0POT«S*C0SZ 

IF ( Z. LT.60. )   GO  TO  10 

IF( Z.LT.70. )   GO  TO  11 

IF  (Z.LT.80. )   GO  TO  12 

IF ( Z. LT.85. )   GO  TO  13 

OP    ■  -567+6. 8*Z 

GO  TO  20 

10  OP   ■  l.+0.0167*Z 
GO   TO  20 

11  OP   »  -*.+0.1*Z 
GO  TO  20 

12  OP   «  -18.+0.3+Z 
GO   TO  20 

13  DP   »  -7*.+Z 

20  0  ■QP0T*(1./(2.71628**(A8SK*0P) ) ) 

IF ( Z. LT.30. )  GO  TO  30 
IF ( Z.LT. 50. )  GO  TO  31 
IFIZ.LT.70.)  GO  TO  32 
IFIZ.LT.80.)  GO  TO  33 
DRATIO  «  -5.8+0.08*Z 
GO   TO  *0 
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30  DR  ATI  0  • 
GO  TO  40 

31  DRATIO  a 
GO   TO  AO 

32  DR ATI  0  - 
GO    TO  40 

33  DR  ATI  0  » 
D  ■ 


1,0 


0.12+0.000667*Z 


0. 065*0. 0025*Z 


-0. 135+0. 0065*Z 


50 


-1.64+0.028*Z 
Q*DRATIO 
IF ( SLOPE. FQ.O. I   GO  TO  50 


'<  (  PI-SLOPE  )  /PI  ) 


QP0T»S*(X4*C0S(H+X5)+X6) 
IFIQPOT.GT.O.)    GO  TO  50 
QPOT   «  0. 


0  *QP0T*(1./(2.716  28**(ABSK*0P>>) 

QTOTAL   •    Q  +  D 


SA    -   90. -Z 
61   SUM1   -    SUM1   ♦  OPOT*DELT 
SUM2   •  SUM2   ♦  0*OELT 
SUM3   •  SUM3  ♦  D*D EL T 

100  CONTINUE 

101  CONTINUE 

SUM*,   •  SUM2   ♦  SU"3 
103  CONTINUE 

CALL  CLOUD 

RETURN 

END 


Subroutine  CLOUD 


C  SUBROUTINE    ACCOUNTS    FOP    PERCENT   CLOUD   COVFR    IN   CONVERTING  POTENTIAL 

C  SOLAR   RADIATION   TO   A  NET  VALUE. 

COMMON  AC TD ATE, ACTUAL! 21 ) »AVETEMC,AVETEMF»BASTENF,ENGBAL»HOLDCAP 
COMMON  CALAIR, CALDEF, CALORIE»CALSNOW,COVDEN, DATE, DATES(3)»DEN 
COMMON  DENSITY,IFIRST, F0OTN0T(16),FREEWAT.IDATE(3  72),ISNOW»ITABLE 
COMMON  KOUNT, LAST USD, LINES, OBSWEOV, RADIN, RADLWN, PADS WN, REFLECT 
COMMON   PACKTEM, PARTICf,PASTINT,PL0T08S,°L0TWE>PRECIP»PREWEaV 
COMMON  SNOMELT,SOBSEOV(372),SPRECIP<3  72),SPREOV(3  72>,SU3TITL(8) 
COMMON  TOTPREC, TITLE (F ) , THR SHL D . TE MP M I N , TE MPMAX , TCOE F F , US E M E AN 
COMMON  X  MAX, USE  POT, SLOPE, A S PECT , XL AT , DELH , ABSK , CL COV , Z AT. SUM  4 
COMMON  KKPP.XLATE 

INTEGER   AC TDATt , DATE, DATES, FOOTNOT, P AS T IN T , PL OTOB S , PLOTWE. SUB TI TL 
INTEGER   TITLE, US E ME AN, US E POT 
DIMENSION  C0LATITI7) 


DICTIONARY 
CLCOV  -  AVERAGE 


CLOUD  COVEP   FOR   THE  DAY 


BASED  ON  LATITUDE. 
38,  .41,  .42,  .44, .52/ 


COLAT   -  A  COEFFICIENT 

DATA  C0LATIT/.35..34,. 
40  ZAT  •  XLATE/10 

DO  10  M«2,7 

IF ( ZAT.LE.M)   GO  TO  20 
10  CONTINUE 
20  N    *  M-l 

COLAT   ■    < COLATIT( D-CPLAT IT (N >  I ♦ < ZAT-N) +COLATIT  (N ) 
30   RADIN   *   SUM4*(1.0-(  (  1-C  OL  A  T  )  *C  L  C  OV  )  ) 
RE  TURN 
END 


Subroutine  DIFTEM 


C  THIS   SUBROUTINE    IS   USED    IN   SIMULATING   CHANGES    IN   SNOWPACK   TEM  — 

C  PERATURES   DURING  PERIODS   WHERE   A  CALORIE  DEFICIT  EXISTS. 

C  DICTIONARY 

C  DEN  -  DENSITY  OF  SNOWPACK 

C  DEPTH  -   DEPTH   OF  SNOWPACK 

C  DIFUS  -  DIFFUSIVITY  OF   THE  SNOWPACK 

C  TAU  -  TIME 

C  OTAU  -  TIME  INTERVAL 

C  FF(I)   -  AN   ARRAY   USED   IN   SIMULATING   SNOWPADK   TEMPERATURES    AT  DIF- 

C  FERENT  DEPTHS. 

C  TEMPK  -  NORMALIZED   M I M  MUM   AIR  TEMPERATURE 

C  T(I)   -  TEMPERATURES   WITHIN  THE  SNOWPACK 

COMMON  ACTDATE,ACTUAL(21),AVETEMC,AVETEMF,8ASTEMF,ENGBAL,HOLDCAP 
COMMON  CALA IR,CALDEF, CALORIE, CALSNOW,COVDEN, DATE, DATES! 3 ) , DEN 
COMMON  DENSITY, I  FIRST, FOOTNOT! 16) , F R E EW AT . I D A TE (  3 72 ) , I SN OW , IT  A B L E 
COMMON  KOUNT, LASTUSD, LINES, OB SWEOV, RADIN, RADLWN,RADSWN, REFLECT 
.  COMMON  PACKTEM, PARTICE,  P AS T I N T, P LOTOBS , PLOTWE, PRECIP, PPEWEQV 
COMMON  S NOME L T» S0BSEQV(372)»SPRECIP!372),SPRE0V(3  72),SUBTITL(8) 
COMMON  TOTPREC,TITLE(fi),THRSHLD,TEMPMIN,TEMPMAX,TCOEFF,USEMEAN 
COMMON  X  MAX, USE  POT, SLOPE  ,  ASPECT, XL  AT, DELH, ABSK, CLCOV, ZAT, SUM*. 
COMMON  KKPP,XLATE 

INTEGER  AC TDATE, DATE, C ATE S. FOOTNOT, PAS TI NT, PL OTOB S, PLOTWE, SUBTITL 
INTEGER  TITLE, US E ME  AN, US E POT 

DIMENSION  DENN(4),F(50)  ,  8  ET  A  (  50  )  ,  GA  M  MA  ( 5  0  )  ,  FF  (  50 )  ,T(  50),  0(50) 
DATA   DENN/.25, .  30, .35,  .40/ 
C  CHECK  TO  SEE   IF   DIFTEK   IS   A  CONTINUING  ESTIMATE  OF  TEMPERATURE  OF 

C  IF   IT  HAS   BEEN  REINITIALIZED. 

210   IF(IFIRST)  10,20,10 
C  CALCULATE  SNOWPACK  DENSITY,   BASED  ON  TIME   OF  YEAR. 

20  IF(DATES(1).LT.5)    GO  TO  30 
IF(DATES(1).GT.9>   GO  TO  40 
DEN*  .40 
GO  TO  50 
40  DEN-  .20 
GO   TO  50 
30  DEN-DENN! DATES ( 1 ) ) 
C  CALCULATE   SNOWPACK  DEPTH   FROM  DENSITY 

50  DEPTH*(PREWEQV/DEN)*2.54 
C  CALCULATE  DIFFUSIVITY 

DIFUS»(0.01)/((2.751-CEN)*0.48) 
C  CALCULATE  TAU  AND  DTAU 

TAU»DIFUS*8 6400 /DEPTH** 2 
DTAU»TAU/48 

C  DETERMINE   THE   AVERAGE  PACK   TEMPERATURE   FROM  THE  EXISTING  CALORIE 

C  DEFICIT 

PACTEM1* (-CALDEF/ (PREWEQV*1. 27 >) +273 
C  SET  THE   INITIAL  TEMPERATURE   ARRAY  AT  THE  AVERAGE  PACK  TEMPERATURE 

DO  60  I«l,50 
60  FF  (  D-PACTEM1 

IFIRST-1 
10  RATIO  •DTAU/(48.0**2) 
A--RATIO 


B*1.0+2.0*RATI0 
C  CONVERT   MINIMUM  AIR   TEMPERATURE   TO  DEGREES  KEiVIN 

TEMPK»((AVETEMF-32.0)*.55555)+2  73 
C  NORMALIZE   EXISTING   PACK  TEMPERATURES 

IFITEMPK. IE. 274.0. AND. TEMPK. GE. 272.0)    TEMPK  «  272.0 

3  DO   7  1*1,50 

7  F( I  )*(273— FF( I ) )/(273-TEMPK) 
DO    4  I»l,49 

4  Til  >*F(I) 
TAU'0.0 

C  PERFORM  CALCULATIONS  OVER   SUCCESIVE  TIME-STEPS 

5  TAU-TAU+DTAU 

IF ( TAU.GT .48*DTAU )    GO  TO  100 
C  SET  BOUNDARY  VALUES 

T ( 2  )  »0.0 
T( 49) «1 .0 

C  COMPUTE   RIGHT-HAND  SIDE'  VECTOR  V 

DO   15  1*2,48 
15  D( I )*T( I) 

D(2)*D(2)+RATI0*T(1> 

D(48)«D(48>+RATI0*T<49) 
C  COMPUTE   NEW  TEMPERATURES 

C  COMPUTE   INTERMEDATE   ARRAYS   BETA   AND  GAMMA 

BETA(2)*B(2) 
GAMMA (2)"D(2)/BETA(2 ) 
DO  6  1*3,48 

BETA(I)»B-A**2/BETA(I-1) 

6  GAMMA ( I ) • (D( I )-A*GAMMA( 1-1 ) ) / BETA! I ) 
C  COMPUTE   FINAL   SOLUTION  VECTOR  T 

T( 48) -GAMMA (48 ) 
DO   11  K»l,46 
I-48-K 

11  T(I)=GAMMA(I)-A*T(I+1)/BETA(I) 
DO   8  1*1,49 

8  F  (  I  )  »  T  (  I  ) 
GO   TO  5 

100  DO   110  1*1,49 

110  FF (  I) *-F( I )*(  273-TEMPK )  +273 

TOT-0.0  TU  _T 

C  ESTABLISH  AN  ARRAY  OF  FINAL   PACK   TEMPERATURES   TO   BE   USED  IN  51 

C  THE  CALCULATIONS   OF   PACK   TEMPERATURES  FOR   THE  NEXT  DAY 

DO   80   1*1, MP1 

80  TOT-TOT  +  FF  (  I  )  -273  tnlT|,.f 
C  DETERMINE   THE  CALORIE  DEFICIT  AND  NEW  AVERAGE  PACK  TEMPERA TU't 

AV-T0T/MP1 

CALDEF»-AV*( PREWEQV*1.27) 

PACKT  EM»A V 

IF (CALDEF )  90,90,95 
90   PACKTEM    "  0.0 

CALDEF*  0.0 
95  RETURN 

END 
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